% the begining of span must be more than P_c

function [p_m_vec, coef_m_mat] = pmvec_coefmmat_of_eigenmodes_analytically_other_ver(n, m, CONSTS)

    k0 = CONSTS.k0;
    a = CONSTS.a;
    eps = CONSTS.eps;
    eta = CONSTS.eta;

    bessel_zeros = pi*(2*n+m+1/2)/2;
    p_m_vec = bessel_zeros*sqrt(eps/abs(eta))*(1/(k0*a));
    coef_m_mat = zeros(4, size(p_m_vec));

    for i = 1:size(p_m_vec, 1)

        p_i = p_m_vec(i); 
        mat = matrix_of_boundary_conditions(p_i, m, CONSTS);
        det(mat);
        %rank(mat);
        K = mat(1:3, 1:3);
        det(K);
        %rank(K);
        C2 = det(K); %(1/1i)*det(K); %det(K);
        f = - C2 * mat(1:3, 4);
        u = K^(-1) * f;
        B1 = u(1);
        B2 = u(2);
        C1 = u(3);

        coef_m_full = [B1; B2; C1; C2];

        % clear small value

        coef_m = coef_m_full; %clear_small_values(coef_m_full, 1e-5);

        eq = mat * coef_m;

        if(false)
            fprintf('equations for p_mn  %f \n', eq);
        end

        epsilon = 1e-3;
        check = 0;
        for ii = 1:size(eq, 1)
            if(abs(eq(ii))<=epsilon)
                check = check + 1;
            end
        end

        %if (check == 4)
            p_m_vec(i) = p_i;
            coef_m_mat(:,i) = coef_m;
        %end

    end
    
end